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Abstract 

Background: Orthopteran migratory locust, Locusta migratoria, and lepidopteran Asian corn borer, Ostrinia furnacalis, are 
two types of insects undergoing incomplete and complete metamorphosis, respectively. Identification of candidate genes 
regulating wing development in these two insects would provide insights into the further study about the molecular 
mechanisms controlling metamorphosis development. We have sequenced the transcriptome of 0. furnacalis larvae 
previously. Here we sequenced and characterized the transcriptome of L migratoria wing discs with special emphasis on 
wing development-related signaling pathways. 

Methodology/Principal Findings: lllumina Hiseq2000 was used to sequence 8.38 Gb of the transcriptome from dissected 
nymphal wing discs. De novo assembly generated 91 ,907 unigenes with mean length of 61 0 nt. All unigenes were searched 
against five databases including Nt, Nr, Swiss-Prot, COG, and KEGG for annotations using blastn or blastx algorithm with an 
cut-off E-value of 10~ 5 . A total of 23,359 (25.4%) unigenes have homologs within at least one database. Based on sequence 
similarity to homologs known to regulate Drosophila melanogaster wing development, we identified 50 and 46 potential 
wing development-related unigenes from L migratoria and O. furnacalis transcriptome, respectively. The identified 
unigenes encode putative orthologs for nearly all components of the Hedgehog (Hh), Decapentaplegic (Dpp), Notch (N), 
and Wingless (Wg) signaling pathways, which are essential for growth and pattern formation during wing development. We 
investigated the expression profiles of the component genes involved in these signaling pathways in forewings and hind 
wings of L migratoria and 0. furnacalis. The results revealed the tested genes had different expression patterns in two 
insects. 

Conclusions/Significance: This study provides the comprehensive sequence resource of the wing development-related 
signaling pathways of L migratoria. The obtained data gives an insight into better understanding the molecular 
mechanisms involved in the wing development in L. migratoria and O. furnacalis, two insect species with different 
metamorphosis types. 
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Introduction 

Insects are the only group of invertebrates that have evolved 
flight [1]. Their wings serve not only as organs of flight, but also 
may be adapted variously as protective covers [2], thermal 
collectors [3], gyroscopic stabilizers [4], sound producers [5], or 
visual cues for species recognition and sexual contact [6] . For those 
insects with wings during the adult stage, complete wings are not 
always visible throughout the life cycle. The insects undergoing 
incomplete metamorphosis have represented functional wings 
during the stage of nymph [7], while the insects going through 
complete metamorphosis only have wing discs inside the body in 



the larval stage [7] . Therefore, the comparative studies on how the 
insect wings are developed will be helpful to understand the insect 
metamorphosis. Additionally, given the high evolutionary conser- 
vation of the proteins involved in the wing development, the 
understanding of the molecular mechanisms involved in wing 
development also sheds light on organogenesis, tissue homeostasis, 
human disease and so on [8-11]. 

Insect wing development is controlled with amazing precision 
and complication. Current understanding of insect wing develop- 
ment mechanisms is mainly from the fruit fly Drosophila 
melanogaster [12-14]. The adult wing of fruit fly is derived from 
the wing imaginal disc, formed at the end of embryonic 
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Figure 1. Observation and comparison of dissected wing discs and adult wings from L. migratoria and O. furnacalis. Note that 0. 

furnacalis adult wings were manually removed scales. 

doi:10.1371/journal.pone.0106770.g001 



development [15]. In Drosophila, the wing imaginal disc is 
subdivided into anterior (A) and posterior (P) compartments at 
very early stage and then further subdivided into dorsal (D) and 
ventral (V) compartments at second instar stage [16-18]. 
Organizers located in the A/P and D/V boundaries coordinate 
the patterning of the wing disc by secreting signal molecules 
including the long-range morphogens Decapentaplegic (Dpp; the 
vertebrate homolog of which is TGFfi) and Wingless (Wg; the 
vertebrate homolog of which is Wnt) [11,19,20], and short-range 
morphogen Hedgehog (Hh) [21]. The produced morphogens form 
gradients to regulate the expression of target genes and control all 
aspects of wing development at cell level via the specific signaling 
pathways [19,22-27]. These signaling pathways are highly 
evolutionarily conserved in many different animals. Some key 
genes for the wing disc development have also been identified in a 
few limited insect species, such as Tribolium castaneum [28]. 
However, knowledge about the identification and involvement of 
wing development-related genes in various insects, especially in 
non-model insects, is still unclear and incomplete. 

Orthopteran migratory locust, Locusta migratoria manilensis, 
and lepidopteran Asian corn borer, Ostrinia furnacalis (Guenee), 
are two different types of insects undergoing incomplete and 
complete metamorphosis, respectively [29,30], Identification of 
candidate genes regulating wing development would provide 
insights into the further study about the molecular mechanisms 
controlling metamorphosis development. Traditionally, such gene 
identification in "non-model" insects relied on degenerate PCR, 
which is labor-intensive, expensive, prone to failure, and only 
produces incomplete fragments [31]. The introduction of novel 
high throughput sequencing technologies greatly facilitates the 
global analysis of the blueprint of development-related genes [32]. 
This technology has been used, for example, to characterize the 
wing development-related genes in the milkweed bug Oncopeltus 
fasciatus [33], the oriental fruit fly Baclrocera dorsalis [34], and 
the salt marsh beetle Pogonus chalceus [35] etc. 



In this study, we combined the Illumina sequencing and de novo 
assembly to obtain and characterize the transcriptome of the wing 
discs of L. migratoria nymph. 91,907 unigenes were assembled 
and 23,359 ones were annotated to known databases. We also re- 
characterized the previous transcriptome of 0. furnacalis larvae 
[36], with special emphasis on wing development-related genes. 
Overall, we identified 50 potential wing development-related 
unigenes from L. migratoria transcriptome, and 46 unigenes from 
O. furnacalis transcriptome, respectively. Additionally, we per- 
formed qRT-PCR analysis to investigate the gene expression 
profiles of several key wing development genes during the stage of 
rapid growth in L. migratoria and O. furnacalis. All these results 
provide valuable information for studying the molecular mecha- 
nisms involved in the insect wing development, and are useful 
resources for further exploring the mechanism how signaling 
pathways control wing development. 

Materials and Methods 

Insect rearing 

Migratory locust L. migratoria manilensis was reared on fresh 
wheat seedlings at 28-30°C, 60% relative humidity with an 18L/ 
6D photoperiod. Asian corn borer (O. furnacalis (Guenee)) larvae 
were reared on an artificial diet at 28°C under a relative humidity 
of 70-90% and a photoperiod of 18L/6D [36]. 

Dissection and total RNA extraction of wing discs 

The wing discs of L. migratoria fifth instar nymph were 
dissected along the wing root with small scissors under microscope. 
Thirty pairs of dissected wing discs were combined, and total RNA 
samples were prepared using TRizol Reagent (TIANGEN, 
Beijing, China) following the manufacturer's instructions. Total 
RNA was dissolved in H 2 0 and RNA quantity was determined on 
a Nanodrop ND-2000 spectrophotometer (NanoDrop products, 
Wilmington, DE, USA). RNA integrity was checked on Agilent 
2100 BioAnalyzer (Agilent Technologies, Englewood, CO, USA). 
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Figure 2. Homology analysis of L. migratoria unigenes. (A) E-value distribution. (B) Similarity distribution. (C) Species distribution. All unigenes 
that had BLASTX annotations within the NCBI Nr database with a cut-off E-value of 10-5 were analyzed. The first hit of each sequence was used for 
analysis. 

doi:1 0.1 371 /journal.pone.01 06770.g002 



Library construction and lllumina sequencing 

Ten (o.g of total RNA was used to isolate mRNA using oligo(dT) 
magnetic beads. The cDNA library was constructed using 
NEBNext mRNA Library Prep Reagent Set (NEB, Ipswich, 
MA, USA) following the manufacturer's protocols. Briefly, 
enriched poly(A) RNA of each sample was fragmented into 200- 
700 nt pieces with RNA Fragmentation Reagents. The cleaved 
RNA fragments were transcribed into the first-strand cDNA using 
random hexamer-primers, followed by second-strand cDNA 
synthesis. The resulting double-stranded cDNA (dsDNA) was 
purified with QiaQuick PCR extraction kit (Qiagen, Hilden, 
Germany) and dissolved in EB buffer. The purified dsDNA was 
treated with T4 DNA Polymerase and T4 Polynucleotide Kinase 
for end-repairing and dA-tailing. After that, they were ligated to 
sequencing adaptors with barcode using T4 DNA ligase. Finally, 
fragments with around 200 bp-length were purified with Qia- 
Quick GelPurify Kit (Qiagen, Hilden, Germany), and used as 
templates for PCR amplification to create the cDNA library. The 
library was sequenced on lllumina HiSeq 2000 (lllumina, San 
Diego, CA, USA) in Baimaike company (Beijing, China). 

Assembly and annotation of transcriptomes 

Raw reads were filtered to remove low quality reads with Q20 
less than 20 and the sequence reads containing adapters and poly- 
A/T tails. The resulting clean reads were assembled to produce 
unigenes using the short reads assembling program - Trinity [37]. 
For functional annotations, we first searched all unigene sequences 
against various protein databases such as Nr, Swiss-Prot, COG, 
and KEGG using BLASTX, and then searched nucleotide 
database Nt using BLASTN, with an E-value cut-off of 10 - ' 



[38]. The BLAST results were used to extract coding region 
sequences (CDS) from the unigene sequences, and translate them 
into peptide sequences. When a unigene happened to have no 
BLAST hits, ESTScan software [39] would be used to determine 
the sequence direction. In addition, we performed the Gene 
Ontology (GO) annotations for each unigene with Blast2GO 
program according to the GO association done by a BLASTX 
against the Nr database [40,41]. 

Identification and sequence analysis of wing 
development-related genes from L migratoria and O. 
furnacalis transcriptome 

The available wing development-related gene sequences from 
Drosophila were used as references to screen L. migratoria 
transcriptome database obtained above and 0. furnacalis 
transcriptome obtained previously [36]. The potential candidates 
of L. migratoria and 0. furnacalis wing development-related genes 
were confirmed by searching the BLASTX algorithm against the 
non-redundant (nr) NCBI nucleotide database using a cut-off E- 
value of 10~ 5 . 

For the sequence analysis of putative wing development-related 
genes identified above, the deduced protein domains were 
determined by using Pfam (http://www.sanger.ac.uk/Software/ 
Pfam/) and SMART (http://smart.embl.de/). Analysis of de- 
duced amino acid sequences, including prediction of signal 
peptide, molecular weight and isoelectric point, was carried out 
in the EXPASY (Expert Protein Analysis System) proteomics 
server (http://www.expasy.org). Sequence comparisons and phy- 
logenetic analysis were performed by MEGA5 software [42]. 
Phylogenetic trees were constructed by the neighbor-joining 
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Figure 3. Gene ontology (GO) assignment for the L. migratoria 
transcriptome. GO assignments (level 2) as predicted for their 
involvement in (A) biological processes, (B) cellular components, and 
(C) molecular functions. The number of unigenes assigned to each GO 
term is shown behind semicolon. 
doi:1 0.1 371 /journal.pone.01 06770.g003 

method, with statistical analysis by the bootstrap method, using 
1000 repetitions. 

Expression assay of several identified wing development- 
related genes 

To investigate the expression profiles of several key wing 
development-related genes in forewings and hind wings, we 
dissected the wing discs as described above, and extracted total 
RNA independently from 3 biological replicates. DNase I-treated 



RNA (lug) was converted into first-strand cDNA using TIAN- 
Script RT Kit (TIANGEN, Beijing, China). The cDNA products 
were diluted 2 fold for use as template. Specific primers for each 
gene were designed and listed in Table SI. L. migratoria Actin and 
O. furnacalis ribosomal protein L8 [rpL8) was used as an internal 
standard to adjust the template amounts in a preliminary PCR 
experiment. The qRT-PCR was performed on a Applied 
Biosystems 7500 Real-time PCR system (Life Technologies, 
Grand Island, NY, USA) using the GoTaq qPCR Master 
(Promega, Madison, WI, USA), according to the manufacture's 
instructions. The thermal cycling conditions for qRT-PCR and 
calculation methods were same as described previously [36]. 

Results and Discussion 

Dissection and observation of wing discs in L migratoria 
and O. furnacalis 

All insects in the Pterygota undergo metamorphosis from 
immature to adult. Among them, insects with incomplete 
metamorphosis, such as migratory locust L. migratoria, have 
young nymph resembling the adult with visible forewings and hind 
wings. Meanwhile, insects with complete metamorphosis, such as 
corn borer O. furnacalis, go through four stage processes from 
egg, larva, pupa, and adult in which only adult has visible wings. 
Although there are some minor differences, the overall appear- 
ances of forewings and hind wings in L. migratoria nymph and 
adult are similar (Fig. 1). However, actual wings are only visible in 
O. furnacalis adult. In fifth-instar O. furnacalis larva, a pair of 
pea-like wing discs composed of large amount of cells was dissected 
from the position where the forewing and hind wing will be 
derived. There is no similarity on the appearance between the 
wing discs and adult wings (Fig. 1). The huge difference about the 
wing development in L. migratoria and O. furnacalis suggests that 
the molecular mechanisms controlling the wing development in 
insects with incomplete or complete metamorphosis might be 
largely different. As a first step to understand these molecular 
mechanisms, it is important and necessary to identify as many as 
possible genes functioning in the wing development. Previously, we 
have combined the Illumina sequencing and de novo assembly to 
obtain the high quality transcriptome from O. furnacalis larvae 
[36]. In this study, we also obtained the data for L. migratoria 
transcriptome, and we can comparatively characterize both 
transcriptomes and identify wing development-related genes. This 
work will provide useful information for studying the molecular 
basis involved in the wing development in L. migratoria and O. 
furnacalis, two insect species with different type of metamorphosis. 

Sequencing and unigene assembly of L migratoria 
transcriptome 

In order to obtain detailed information about L. migratoria 
transcriptome, we prepared cDNA from the wing discs of fifth- 
instar nymph, and subjected it to Hiseq 2000 sequencing. After 
cleaning of dirty reads and quality checks, a total of 83,000,540 
high-quality clean reads (SRA accession number SRX491784) 
with a cumulative length of 8,382,624,748 nucleotides (8.38 Gb) 
were generated from L. migratoria wing disc library. The GC 
percentage of the reads is 42.65%, which is comparable with 
genome sequence of other insects. Using Trinity software by the 
manner of paired-end joining and gap-filling, these reads were 
assembled into 91,907 unigenes longer than 200 nt (mean length 
of 610 nt and N50 of 1024 nt). The size distribution indicated that 
the ratio of unigenes with a length of 200-1000 bp was 86.49%, 
while the length of 12,420 (13.51%) unigenes was more then 
1000 bp (Fig.Sl). It was significantly larger than that in previous 
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Figure 4. Clusters of orthologous groups (COG) classification of L. migratoria unigenes. A total of 6,508 produced functional annotations 
were among the 25 categories. The Y-axis shows the number of unigene in each COG term. 
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insect transcriptome projects [43,44]. The assembled sequences 
have been deposited in the NCBI Transcriptome Shotgun 
Assembly (TSA) Database under the accession GBDZ00000000. 

Functional annotation and classification 

In order to annotate the unigenes, the obtained sequences were 
first aligned by BLASTX to various protein databases of nr, Swiss- 
Prot, KEGG, COG, and GO (E-value<10" 3 ), and then aligned 
by BLASTN to nucleotide database nt (E-value<10 5 ). Among 
91,907 unigenes, 20,746 (22.6%), 8,737 (9.5%), 11,450 (12.5%), 
5,312 (5.8%), 4,986 (5.4%), and 10,644 (11.6%) ones were 
annotated in nr, nt, Swiss-Prot, KEGG, COG, and GO, 
respectively (Table S2). A total of 23,359 unigenes (25.4%) were 
annotated to at least one database. The remaining 68,548 
unigenes (74.6%) were not annotated to any referred databases. 
The low annotated percentage might be due to the transcripts 
derived from the cDNA of untranslated regions, assemblage 
errors, and nonconserved areas of proteins where homology is not 
detected [43]. Another possibility was that a large part of the genes 
in L. migratoria transcriptome database were with unknown 
functions, and the un-annotated unigenes were the potential 
sources of novel genes. 

The functional annotations of unigenes were performed mainly 
based on the BLASTX results against the nr database. Among the 
20,746 unigenes annotated to nr database, 7,883 (38%) showed 
strong homology (E-value smaller than le-50) (Fig. 2A). The 
identity comparison showed 5,809 (28%) unigenes have more than 
60% identity with other insects (Fig. 2B). For species distribution 
of the top BLAST hits against the Nr databse, L. migratoria 
unigenes had the greatest number of matches with Tribolium 
castaneum (2,490 unigenes, 12%), followed by Pediculus humanus 
corporis (1,452 unigenes, 7%), and Acyrlhosiphon pisum (1,245 
unigenes, 6%) (Fig. 2C). 

Gene Ontology (GO) assignment programs were used to classify 
the functions of all unigenes. GO contains three categories: 
biological process, cellular component, and molecular function. 



Using the Blast2GO and WEGO software, 43,415, 24,188, and 
13,706 unigenes were associated with biological process (22 sub- 
categories), cellular component (16 sub-categories), and molecular 
function (17 sub-categories), respectively (Fig. 3). Among the 
biological process assignments, cellular processes (15.0%) and 
metabolic processes (14.4%) represented the most abundant 
subcategories. It indicated the importance of cell cycle, generation 
as well as metabolic activities in wing development stage. Under 
the category of cellular component, the top 3 sub-categories were 
cell part (21.0%), cell (20.2%) and organelle (16.0%). In the 
molecular function category, binding (43.0%) and catalytic 
activities (36.9%) were the most abundant (Fig. 3). The subcate- 
gories taking up the largest two proportions in each category were 
consistent with that in transcriptomic studies of other insects 
[45,46]. 

The unigenes were further compared against COG for the 
analysis of the putative protein functions. A total of 6,508 unigenes 
were functionally classified into 25 COG categories (Fig. 4). The 
cluster for 'General function prediction only' (1,407 unigenes) 
constituted the largest group, followed by 'Replication, recombi- 
nation and repair' (809 unigenes), 'Translation, ribosomal 
structure and biogenesis' (560 unigenes), and 'Transcription' 
(441 unigenes) (Fig. 4). The group of 'signal transduction 
mechanisms' contains 354 unigenes, the fifth largest one within 
the COG classification of the L. migratoria transcriptome. It 
suggests the possible importance of signaling transduction path- 
ways in wing development in L. migratoria. 

Identification and comparison of wing development- 
related genes in L. migratoria and O. furnacalis 
transcriptome 

Wing development is a complicated biological process compris- 
ing multiple integrated signaling pathways. In particular, four 
signaling pathways - Notch, Hedgehog (Hh), Decapentaplegic 
(Dpp), Wingless (Wg) - play important roles in growth control and 
cell fate determination during the wing development 
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[20,22,25,47-49]. In order to obtain the global perspective on the 
wing development-related transcripts in L. migratoria, we 
searched the assembled transcriptome for orthologous genes 
known to be involved in the above four pathways in Drosophila. 
In total, we have identified 50 unigenes with significant similarity 
to Drosophila wing development-related genes, including 7 ones in 
Notch pathway, 12 in Hh pathway, 9 in Dpp pathway, and 10 in 
Wg pathway (Table 1). Meanwhile we performed the similar 
analysis for 0. furnacalis transcriptome reported previously [36], 
and identified 46 unigenes with high potential to be related to the 
wing development (Table 1). 

Genes involved in the Hh signaling pathway. The Hh 
pathway is an important signaling cascade in Drosophila to pattern 
the embryonic cuticles and adult appendages [50]. It is also vital 
for diverse aspects of animal development and essential in humans 
to pattern their limbs and internal organs [51,52]. The pathway 
takes its name from its polypeptide ligand, an intercellular 
signaling molecule called as Hh [53]. Hh was initially discovered 
as a segment polarity gene in Drosophila [54]. Mammals have 
three Hh orthologues, Sonic (Shh), Indian (Ihh), and Desert (Dhh) 
[55]. We have identified the potential hh genes, a37032 from L. 
migratoria and Unigene 16362 from 0. furnacalis transcriptome, 
respectively (Table 1). Although these two unigenes were incom- 
plete, they have 66% and 73% amino acid sequence identity to 
Drosophila Hh (Fig. 5). On the other hand, Drosophila Hh is 
synthesized as inactive precursors with an N-terminal signaling 
region linked to a C-terminal autoprocessing region which begins 
with a cysteine for an acyl rearrangement analogous to step 1 of 
the protein splicing pathway [56]. This Cys residue is completely 
conserved in both L. migratoria and O. furnacalis Hh homologues 
(Fig. 5). The conserved motif A, B, F, J, K, and L in Drosophila 
Hh [57] were found in L. migratoria, but not in 0. furnacalis Hh 
because its C-terminus was presently unknown (Fig. 5). The 
phylogenetic analysis of Hh from different insect species and other 
animals reveals that almost all nematode Hh-related proteins form 
a distinct clade while the other eukaryote Hh proteins, including 
L. migratoria a37032 and O. furnacalis Unigene 16362, cluster 
into another Hh clade (Fig. 6). L. migratoria a37032 is grouped 
with other Orthopteran Hhs, and 0. furnacalis Unigene 16362 is 
clustered with other Lepidopteran Hhs. We performed qRT-PCR 
analysis to analyze the transcriptional expression of hh during 
different developmental stages. The mRNA abundance of hh in L. 
migratoria remained unchanged from the first instar through the 
fifth instar stage. The transcript abundance of 0 . furnacalis hh was 
consistent before the fourth instar stage, but increased significantly 
in the fifth instar larvae (Fig. 7). This result also suggests that the 
quality of the transcriptome assembly is high enough to be used for 
primer design. 

Hh signaling is mediated by a multi-component receptor 
complex in the cell membrane. This receptor complex consists 
of a 1 2-span transmembrane protein, Patched (Ptc) as the receptor 
and a 7-span transmembrane protein, Smoothened (Smo) as the 
obligatory signal transducer across the plasma membrane [51,58]. 
When extracellular Hh binds to and is inhibited by Ptc, Smo starts 
to accumulate, and inhibit the proteolytic cleavage of zinc-finger 
transcription factor Cubitus interruptus (Ci) which is normally 
bound by the kinesin-like protein Costal-2 (Cos2). The intact Ci 
protein then translocate into the nucleus, allowing the transcrip- 
tion of some genes such as dpp. In the absence of Hh, Ptc blocks 
Smo activity, and full-length Ci protein is degraded into Ci 
fragment (CiR) that functions as a transcriptional repressor to 
block the transcription of target genes (Fig. S3A) [50]. By 
searching L. migratoria and O. furnacalis transcriptome using 
Drosophila corresponding genes as reference, we have identified 
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Figure 5. Alignments of potential Hhs in L. migratoria, O. furnacalis and Drosophila. The deduced amino acid sequences of L. migratoria 
a37032 and 0. furnacalis Unigene16362 were compared with Drosophila Hh (Dm_Hh). Completely conserved amino acids are indicated by "*", and 
conservative substitutions by ":" and "." below the sequences. The conserved cysteine residue for an acyl rearrangement analogous to step 1 of the 
protein splicing pathway is boxed. The conserved motifs in Drosophila Hedgehog are underlined. 
doi:1 0.1 371 /journal.pone.01 06770.g005 



the potential 1:1 orthologs of Ptc, Smo, Ci from both 
transcriptome databases (Table 1). It suggests that Hh signaling 
pathway is conserved in both L. migratoria and O. furnacalis. 

Genes involved in Dpp signaling pathway. Dpp, directed 
by the Hh signaling, is a member of the bone morphogenefic 
protein (BMP) growth factor family [59] . Dpp forms a long-range 
dynamic and precise gradients to control cell survival, cell 
morphogenesis, cell proliferation, cell differentiation during the 
wing development [22,25,26,60]. When a cell receives a Dpp 
signal, its heteromeric receptor complex composed of type I 
receptor Thickveins (Tkv) and type II receptor Punt is activated 
[61,62]. Once activated, the receptors are able to phosphorylate 
an intracellular protein called mothers against Dpp (Mad) [63]. 
The phosphorylated Mad then associate with the Medea (Med), 
and the complex translocates to the nucleus where it binds to DNA 
and activates or suppresses the expression of the target genes in 
conjunction with other transcription factors [64]. Within this 
signaling pathway, the Drosophila Hox gene Ultrabithorax (Ubx) 



restricts both the transcription and the mobility of Dpp [65]. 
Genes activated by Dpp signaling pathway include transcription 
factor oplomotor-blind (omb) and two members of Spalt (sal) family 
spalt major (salm), spalt related (salr) and so on (Fig. S3B) [66,67]. 

In this study, we identified almost all transcripts of the above key 
components in Dpp signaling pathway (Table 1). Among them, 
one unigene with sequence similarity to Drosophila Dpp was 
identified from the transcriptome of L. migratoria (a 136595) and 
O. furnacalis (Unigene7326), respectively. The 5'-end of the 
cDNA sequences of both potential Dpp is incomplete. In order to 
establish the orthology of these two fragments, we performed a 
phylogenetic analysis incorporating a selection of BMP growth 
factor members from Drosophila, other arthropods and a variety of 
deuterostome taxa. The resulting phylogenetic tree distinguishes 
two groups of proteins. One group consists of the arthropod dpp 
genes and their deuterostome homologs, the BMP2/4 genes. The 
second group comprises the Drosophila TGF-beta genes screw 
(sew) and glass bottom boat (gbb), and the remaining deuterostome 
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Figure 6. Phylogenetic analysis of Hhs. The amino acid sequences of L. migratoria a37032 and 0. furnacalis Unigene16362, together with 68 
potential Hhs from other organisms were used to build the neighbor joining tree. L. migratoria a37032 and 0. furnacalis Unigene16362 are boxed. 
Nematode Hh-related proteins and other eukaryote Hh proteins are indicted with brackets. The arrows at nodes denote bootstrap value greater than 
700 from 1000 trials. 
doi:1 0.1 371 /journal.pone.01 06770.g006 



BMP genes with TGF-beta homology (Fig. 8). The identified 
fragments, L. migratoria al36595 and O. furnacalis Unigene7326, 
reside in the dpp group with statistical support (reliability 
value = 98 and 99 for al36595 and Unigene7326, respectively). 
Therefore, they are postulated as potential dpp genes. We 
conducted qRT-PCR analysis to investigate the expression profile 
of dpp. L. migratoria dpp was expressed at high level in the first 
instar nymphs, then decreased significantly in the second, third, 
and fourth instar nymphs, but significantly increased in the hind 
wings of the fifth instar nymph. The mRNA level of O. furnacalis 
dpp remained unchanged in all tested stages (Fig. 7). The different 
expression patterns oidpp in two insects suggest that dpp might be 
related to the different metamorphosis. 

In addition, it is interesting that only one mad gene was present 
in Drosophila but two transcripts encoding for potential mad were 
identified from both L. migratoria and 0. furnacalis transcrip- 
tomes (Table 1). Two Mads in L. migratoria (al 1 173 and al2619) 
shared 60% identity while two O. furnacalis Mads (Unigene 16360 
and Unigene6852) had 61% identity in amino acid sequences. A 
similar situation happened to the identification of med. There is 
only one med gene in Drosophila but two transcript fragments were 
identified for both L. migratoria and 0. f urnacalis med, including 
a26380 and al 165 1 in I. migratoria, and unigene4133 and 
unigene 17363 in O. furnacalis, respectively. The encoded amino 



acid sequences of a26380 and unigene4133 were highly similar to 
the N- terminus of Drosophila Med while al 165 1 and uni- 
gene 17363 were highly similar to the C-terminus of Drosophila 
Med (Fig. S2). We doubted that the two unigenes from each 
transcriptome were just two partial fragments within med gene. 
However, no overlapped sequences were observed in a26380 and 
all 65 1 , or in unigene4133 and unigenel7363 (Fig. S2). One 
possible reason might be due to the sequencing errors because the 
predicted overlapping part was just located at the terminus of 
sequenced fragment. Further experiments are required to deter- 
mine the full length of med in these two insect species. 

Genes involved in the Notch signaling pathway. The 
Notch signaling pathway is highly conserved throughout the 
animal kingdom. It regulates cell-fate determination during 
development and maintains adult tissue homeostasis [47]. The 
key component of this signaling pathway is the Notch receptor. 
Notch is a 300 kDa single-pass transmembrane protein, composed 
of a large extracellular domain, a single transmembrane portion, 
and a small intracellular region [68-70]. After binding to its 
ligands, Delta and Serrate (known as Jagged in mammals), inactive 
Notch precursor undergoes two proteolytic cleavage events: the 
first cleavage is catalyzed by ADAM-family metalloproteases; the 
second cleavage is mediated by y-secretase which is an enzyme 
complex containing presenilin, nicastrin, PEN2 and APH1 [71]. 
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Figure 7. qRT-PCR analysis of transcriptional levels of hh, dpp, and Notch in L. migratoria and O. furnacalis. Actin and ribosomal protein L8 
(rpL8) were used as an internal standard to normalize the templates in L migratoria and 0. furnacalis, respectively. RNA was extracted from wing discs 
of forewings (F) or hind wings (H) collected from the first-instar (L1), second-instar (L2), third-instar (L3) : fourth-instar (L4) and fifth-instar (L5) L 
migratoria nymph or 0. furnacalis larvae. The bars represent the mean ± S. D. (n = 3). Bars labeled with different letters are significantly different (one- 
way ANOVA, followed by the Newman-Keuls test, P<0.05). 
doi:1 0.1 371 /journal.pone.01 06770.g007 



The second cleavage liberates the Notch intracellular domain 
(NICD), which then migrates into the nucleus and cooperates with 
the DNA-binding protein CSL (also named as CBF1, Suppressor 
of hairless (Su(H)), LAG-1, or RBP) and its co-activator 
Mastermind (Mam) to promote the transcription of downstream 
target genes, such as wg (Fig. S3C) [72]. Using the corresponding 
components in Drosophila as referred sequences, we have 
identified 9 transcripts for potential components in Notch signaling 
pathway from L. migratoria and O. furnacalis transcriptomes, 
respectively (Table 1). Four out of nine unigenes in L. migratoria 
are complete while only one unigene is complete in 0. furnacalis. 
The possible reason was that the sequencing quality of L. 
migratoria transcription was better than that of O. furnacalis 
transcriptome (8.38 Gb vs. 4.72 Gb). Given the importance of 
Notch in this signaling pathway, we selected it for further analysis. 
Similar to the case in other invertebrates, we only identified one 
transcript for Notch gene, a 1269 from L. migratoria transcrip- 
tome, and Unigene6905 from 0. furnacalis transcriptome 
(Table 1). The canonical Notch protein consists of three repeated 
sequence motifs: the extracellular domains contain 10 — 36 copies 
of an ~40-amino acid epidermal growth factor-like (EGFL) 
sequence motif and 3 copies of an ~40-amino acid lin/Nolch/glp 
(LNG) sequence motif; the intracellular domains contain 6-7 
copies of a cdclO/SWl 6 /ankyrin (CDC/ANK) sequence motif 
flanked by stretches of nonrepetitive sequences [73]. L. migratoria 



a 1269 encodes a 2,484-amino acids full-length protein which has 
63% similarity to Drosophila Notch. L. migratoria Notch contains 
36 EGFL tandem repeats, 3 LNG repeats, and 7 tandem ANK 
repeats (Fig. 9A), suggesting it is a canonical Notch. 0. furnacalis 
Unigene6905 encodes a 780-amino acids polypeptide with 70% 
similarity to Drosophila Notch. Only 20 tandem EGFL repeats are 
predicted in the current identified Unigene6905 fragment 
(Fig. 9A). It is unknown whether O. furnacalis Notch also contain 
classic LNG and ankyrin repeats because the current transcript is 
incomplete. We performed the phylogenetic analysis for L. 
migratoria and O. furnacalis Notch and 30 Notch sequences 
from other species to investigate the evolutionary relationship of 
the Notch protein family. As shown in Fig. 9B, the Notch from 
insects forms a separate group which includes L. migratoria and 0. 
furnacalis Notch, and the Notch from vertebrates is clustered into 
another group. Notch from Ciona, sea urchin, and amphioxus is 
grouped with vertebrate Notch, however, with a low bootstrap 
value of 82. Caenorhabditis elegans Notch is out of any group, 
showing a great differentiation from the other taxa. We analyzed 
the expression profiles of Notch using qRT-PCR methods. As 
shown in Fig. 7, the expression profile of Notch was similar to that 
of Hh in two insects. It kept unchanged from the first instar 
through the fifth instar stage in L. migratoria, while it increased 
significandy in O. furnacalis fifth instar larvae (Fig. 7). 
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Figure 8. Phylogenetic analysis of Dpps. Except for L migratoria a1 36595 and 0. furnacalis Unigene7326 (marked in red), the used amino acid 
sequences of other 33 Dpps are from mouse (Mm), zebrafish (Dr), lancet (Branchiostoma floridae; Bf), acorn worm (Ptychodera flava; Pf), the sea 
urchins Strongylocentrotus purpuratus (Sp) and Lytechinus variegatus (Lv), fruit fly (Dm), flour beetle (Tribolium castaneum; Tc), sawfly (Athalia rosae; 
Ar), buckeye butterfly (Junonia coenia; Jc), the grasshoppers Schistocerca americana (Sa) and S. gregaria (Sg), cricket (Gryllus bimaculatus; Gb), 
Cupiennius salei (Cs), Ceratitis capitata (Cc), Blattella germanica (Bg), Polyrhachis vicina (Pv), Episyrphus balteatus (Eb), Apis mellifera (Am), Pheidole 
morrisi (Pm), Biston betularia (Bb), Bombyx mori (Bm), Culex quinquefasciatus (Cq), Nasonia vitripennis (Nv). The branches specific for arthropod Dpps 
and TGF-beta homologs are shaded in blur and yellow, respectively. For explanation of the arrows see Fig. 6. 
doi:1 0.1 371 /journal, pone.01 06770.g008 



Genes involved in the wg signaling pathway. Wg (the 
vertebrate homolog of which is Wnt) is transcriptionally activated 
by Notch signaling, and Wg signaling pathway plays critical roles 
in axis patterning, cell fate specification, cell proliferation, and cell 
migration etc [74]. During the development of Drosophila wing, 
canonical Wg signaling pathway specifies pattern formation along 
the dorsal/ventral (D/V) axis while Dpp signaling pathway play 
this role along the A/P axis [75]. The ligand Wg is the founding 
member of the Wnt family, and is a secreted lipid-modified 
signaling glycoprotein that has 350-400 amino acids in length 
[76] . Wg is expressed at the D/V boundary and forms a stable and 
long-range gradient by symmetrically diffusing at both sides of the 
boundary. When Wg binds to a receptor complex consisting of the 
seven-transmembrane protein Frizzled (Fz) and the single-pass 
transmembrane protein Arrow (Arr, homologous to murine and 
human low density lipoprotein (LDL) receptor-related protein 5 or 
6 (LRP5/6)), the downstream cytoplasmic protein disheveled (Dsh 
in Drosophila and Dvl in vertebrates) is activated [7 7] . Dsh in turn 
inhibits glycogen synthase kinase (GSK)-3P in the P-catenin 
destruction complex, which mainly consists of Axin, GSK-3P, 
adenomatous polyposis coli (APC) and P-catenin (armadillo (arm) 
in Drosophila). Consequendy, P-catenin accumulates in the 
cytoplasm, and stabilized P-catenin then translocates into the 
nucleus and acts together with the transcription factor Pangolin to 
regulate the transcription of Wg target genes (Fig. S3D) [77]. In 
this study, we identified 10 and 11 unigenes for the known 
components in the Wg signaling pathway from L. migratoria and 
O. furnacalis transcriptome, respectively (Table 1). No putative 
Frizzled 2 <fz2) ortholog was identified in L. migratoria 
transcriptome. A possible reason is that/z2 gene is missing in L. 
migratoria because of the evolutionary event. The other reason 
with higher possibility is that the transcript level oifz2 is low and it 



is not captured in the RNA-seq. We attempted to perform further 
analysis for identified wg gene, a40139 in L. migratoria and 
Unigenel219 in O. furnacalis. However, unigenel219 only 
encodes a 9 7 -amino acid peptide which is too short to have no 
common sequences with other Wg during the alignment. 
Therefore, we failed to conduct phylogenetic analysis to reveal 
the evolutionary relationship of Wg. Additionally, we identified 
another 1 2 and 9 unigenes from L. migratoria and O. furnacalis 
transcriptome, respectively (Table 1), which were potentially 
involved in the wing development, including homologs to 
Drosophila apterous, engrailed, homothorax, teashirt, epidermal 
growth factor receptor, rhomboid, nubbin, pannier, notum, fat, 
four-jointed, dally-like. 

Conclusions 

In summary, we sequenced and characterized the transcriptome 
from the wing discs of L. migratoria nymph. The assembled 
sequence data comprising 91,907 unique transcripts provides a 
comprehensive sequence source for future L. migratoria study. We 
identified a large set of genes relevant to wing development with 
high significance, especially the genes involved in four signaling 
pathways - Notch, Hh, Dpp, and Wg signaling pathways, from L. 
migratoria transcriptome and another O. furnacalis transcriptome 
obtained previously. The explored wing development-related 
genes constitute an integrated picture of the development network, 
which provides the valuable clues for a better understanding of the 
wing development in L. migratoria and O. furnacalis. These 
development repertoire genes appear to be evolutionarily 
conserved to different extent. Functional analyses are necessary 
to verify our predictions. Nevertheless, the framework of 
information presented in this study should help to further 
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Figure 9. (A) Schematic representation of the L. migratoria and O. furnacalis Notch. The domain organization was predicted using the SMART 
program (http://smart.embl.de/). Question mark means the end is incomplete. (B) Phylogenetic analysis of the Notch from insects, hydrobios and 
vertebrates. The unrooted tree was generated using the distance-based method Neighbor-Joing, with Caenorhabditis elegans protein as outgroup. 
The branch lengths reflect evolutionary divergence. The circled bootstrap value indicates that L. migratoria a1269 and 0. furnacalis Unigene6905 
(marked in red) belong to insect Notch group. 
doi:1 0.1 371 /journal.pone.01 06770.g009 



understand the complex molecular mechanisms involved in wing 
development in L. migratoria and 0. furnacalis, two insect species 
with different type of metamorphosis. 

Supporting Information 

Figure SI Assembled unigene length distribution of L. migra- 
toria transcriptome. The jc-axis indicates unigene size and the y- 
axis indicates the number (left) or percentage (right) of unigenes of 
each size. 
(EPS) 

Figure S2 Alignments of potential Meds in L. migratoria, O. 
furnacalis and Drosophila. Note: L. migratoria a26380 (green) and 
O. furnacalis Unigene4133 (purple) match well with the N- 
terminus of Drosophila Med (blue), while L. migratoria al 1 65 1 
(black) and 0. furnacalis Unigene 1 7363 (red) are highly similar to 



the C-terminus of Drosophila Med (blue). The part without 

identity is shaded. 

(EPS) 

Figure S3 Schematic drawing of the Hh (A), Dpp (B), Notch (C), 
and Wg (D) signaling pathways. (A) In the presence of Hh, Hh 
binds to its receptor Ptc and Ptc relieves Smo repression. Smo 
accumulates and is activated by phosphorylation. This promotes 
its association with the complex including Costal 2 (Cos2), Fused 
(Fu), Suppressor of Fu (Su(Fu)) and Cubitus interruptus (Ci). 
Uncleaved Ci is released from the complex and acts as 
transcriptional activator in the nucleus and induces expression of 
target genes. (B) Dpp binds to its receptor complex including the 
type I receptor Tkv and type II receptor Punt which induces the 
phosphorylation of Tkv and in turn phosphorylates the signal 
transducer Mad. Phosphorylated Mad forms complex with Med 
and then translocates to the nucleus and induces expression of 
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target genes. (C) Notch receptor is a transmembrane protein, with 
a large extracellular domain (NECD), a transmembrane domain 
and an intracellular domain (NICD). When its ligands Dl or Ser 
bind to the Notch, two proteolytic cleavage events are activated: 
the first cleavage is catalyzed by ADAM-family TACE metallo- 
proteases; the second cleavage is mediated by y-secretase which is 
an enzyme complex containing presenilin, nicastrin, PEN2 and 
APH1. The second cleavage liberates the NICD, which then 
migrates into the nucleus and cooperates with the DNA-binding 
protein CSL and its co-activator Mam to promote the transcrip- 
tion of downstream target genes. (D) Wg binds to the receptors Fz 
and Arr. The Wg-Arr-Fz complex binds to and activates the 
protein Dsh which then disrups the complex of APC/Axin/Sgg 
and results in an increase of free cytosolic Arm. Arm then 
translocates to the nucleus where it binds to transcription factor 
Pangolin (dTCF) and activates the expression of target genes. 
(EPS) 
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